1 Introducción

El siguiente documento corresponde al primer avance del trabajo de investigación del curso Visualización de datos y análisis estadístico del Posgrado Analista de Datos de la Universidad Cenfotec.

A continuación, se exponen el tema a tratar en esta investigación, los objetivos planteados, asi como una muestra del dataset que se empleara en dicha investigación.

2 Tema de investigación

“Estadisticas y Modelado de la distribución biogeográfica de Pumas y Jaguares de un sector del continente americano”

4 Objetivos

    1. Estadísticas espacio-temporales sobre los pumas y jaguares de un sector del continente americano, en las que se evalua: densidad por zonas y por pais, patrones de distribución dentro del area de estudio y frecuencia de observaciones por mes
    1. Analisis de interpolación de Krigging respecto del número de individuos observados por año
    1. Modelado de la distribución de los pumas y jaguares del área de estudio, según la máxima entropía con la cual se determina la probabilidad de condiciones adecuadas para la especie

5 Librerias

library(dplyr)
library(sp)
library(DT)
library(leaflet)
library(raster)
library(rasterVis)
library(rgdal)
library(ggplot2)
library(plotly)
library(sf)
library(here)
library(rasterVis)
library(devtools)
library(proto)
library(dygraphs)
library(highcharter)
library(echarts4r)

6 Procesos de extracción, transformación y carga de las Bases de datos

En esta sección se realizan una serie de procesos que conllevan al desarrollo de las bases de datos finales.

6.1 1. Base de datos sobre felinos

df <- read.csv("E:/Analista de datos/DAT004/Proyecto Final/Datasets/samples/wildcats_df.csv", stringsAsFactors=F, sep=",", na.strings=c("","NA"), encoding="UTF-8")

colnames(df)
## [1] "X.U.FEFF.common_name" "scientific_name"      "latitude"            
## [4] "longitude"            "created_at"           "updated_at"          
## [7] "place_country_name"

6.1.1 Cambio de nombre de las columnas

df <- df[,c(1:5,7)]

names(df)[1] <- "Nombre_Comun"
names(df)[2] <- "Nombre_Cientifico"
names(df)[3] <- "Latitud"
names(df)[4] <- "Longitud"
names(df)[5] <- "Fecha_Registro"
names(df)[6] <- "Pais"

colnames(df)
## [1] "Nombre_Comun"      "Nombre_Cientifico" "Latitud"          
## [4] "Longitud"          "Fecha_Registro"    "Pais"

6.1.2 Resumen de datos por especie

6.1.2.1 Conteo inicial

df%>%
  dplyr::group_by(Nombre_Comun)%>%
  dplyr::summarise(length(Nombre_Comun))

6.1.2.2 Correccion de datos en la columna de “Nombre_Comun”

ndf <- df%>%
  mutate(Nombre_Comun=ifelse(Nombre_Comun=="Puma de América del Norte", "Puma",Nombre_Comun))%>%
  mutate(Nombre_Comun=ifelse(Nombre_Comun=="Puma de América del Sur", "Puma",Nombre_Comun))

6.1.2.3 Resumen final de datos por especie

ndf%>%
  group_by(Nombre_Comun)%>%
  summarise(length(Nombre_Comun))

6.1.2.4 Limpieza de NA

ndf <- ndf[complete.cases(ndf), ]

6.1.3 Mapa de la distribicion de las especies a nivel mundial

cf <- colorFactor(c("#ffa500", "#13ED3F"), domain=ndf$Nombre_Comun)
m <- leaflet(df) %>% addTiles(group = "OSM (default)") %>%
  addProviderTiles(providers$Stamen.Toner, group = "Toner") %>%
  addProviderTiles(providers$Stamen.TonerLite, group = "Toner Lite") %>%
  setView(-63.54105072, -8.88123188, zoom = 2) %>% 
  addCircles(~Longitud, ~Latitud, popup=ndf$Nombre_Comun, group = ndf$Nombre_Comun ,weight = 3, radius=4, 
              color=cf(ndf$Nombre_Comun), stroke = TRUE, fillOpacity = 0.5)  %>%
  addLegend("bottomright", colors= c("#ffa500", "#13ED3F"), labels=c("Puma", "Jaguar"), title="Especies") %>% 
  addLayersControl(
    baseGroups = c("OSM (default)", "Toner", "Toner Lite"),
    overlayGroups = ndf$Nombre_Comun,
    options = layersControlOptions(collapsed = TRUE)
  )
m

6.2 2. Bases de datos sobre variables ambientales del área de estudio

El area de estudio con la que se pretende trabajar corresponde a Centro y Sudamrica, asi como el Caribe.

Algunas de las coberturas que se utilizan para conocer el habitat de las especies son variables climáticas que se derivan de los datos proporcionados por el Panel Intergubernamental sobre Cambio Climático y se produjeron utilizando Interpolación de basada en lecturas tomadas en estaciones meteorológicas de todo el mundo desde 1961 hasta 1990.

Los coverturas corresponden a:

*cld6190_ann : cobertura de nubes, anual

*dtr6190_ann : rango de temperatura diurna, anual

*frs6190_ann : frecuencia de heladas, anual

*pre6190_ann : precipitacion, anual

*pre6190_I1 : precipitacion, enero

*pre6190_I4 : precipitacion, abril

*pre6190_I7 : precipitacion, julio

*pre6190_I10 : precipitacion, octubre

*tmn6190_ann : temperatura promedio, anual

*tmp6190_ann : temperatura minima, anual

*tmx6190_ann : temperatura maxima, anual

*vap6190_ann : presion de vapor, anual

r_files <- list.files(here::here("Datasets", "coverages"),
                       full.names = T)
r_files
##  [1] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/cld6190_ann.asc"
##  [2] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/dtr6190_ann.asc"
##  [3] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/ecoreg.asc"     
##  [4] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/frs6190_ann.asc"
##  [5] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/h_dem.asc"      
##  [6] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_ann.asc"
##  [7] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_l1.asc" 
##  [8] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_l10.asc"
##  [9] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_l4.asc" 
## [10] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/pre6190_l7.asc" 
## [11] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/tmn6190_ann.asc"
## [12] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/tmp6190_ann.asc"
## [13] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/tmx6190_ann.asc"
## [14] "E:/Analista de datos/DAT004/Proyecto Final/Datasets/coverages/vap6190_ann.asc"
st <- raster::stack(r_files) 

raster::plot(st)

6.2.0.1 Muestra de una de las variables ambientales: Ecoregiones

ecoreg = raster(r_files[3])
raster::plot(ecoreg, main="Ecoregiones del area de estudio")

raster::hist(ecoreg, main="Histograma de las Ecoregiones del area de estudio", col="green")

6.2.0.2 Delimitacion del area de estudio

A continuacion, se carga el un archivo espacial de tipo vectorial correspondiente a los paises que conforman el area de estudio.

area_estudio <- sf::st_read(dsn="E:/Analista de datos/DAT004/Proyecto Final/Datasets/samples/map.geojson", layer="map", quiet=TRUE)

Y procede a crear una visualizacion del mismo.

map.area <- plot_geo(area_estudio, name=~pais, color = I("#89C5DA"))%>% hide_legend()

map.area

6.3 3. Selección de la base de datos de las especies respecto del área de estudio

coords <- data.frame(x=ndf$Longitud, y=ndf$Latitud)
prj <- CRS("+proj=longlat +datum=WGS84 +no_defs")
points <- SpatialPointsDataFrame(coords, ndf, proj4string = prj)
wildcats <- crop(points,ecoreg)
wildcats <- data.frame(wildcats)
wildcats <- wildcats[,1:6]

6.3.1 Mapeo de los ubicaciones registradas de las especies en el área de estudio

cf2 <- colorFactor(c("#ffa500", "#13ED3F"), domain=wildcats$Nombre_Comun)
m2 <- leaflet(wildcats) %>% addTiles(group = "OSM (default)") %>%
  addProviderTiles(providers$Stamen.Toner, group = "Toner") %>%
  addProviderTiles(providers$Stamen.TonerLite, group = "Toner Lite") %>%
  setView(-63.54105072, -8.88123188, zoom = 2) %>% 
  addCircles(~Longitud, ~Latitud, popup=wildcats$Nombre_Comun, group = wildcats$Nombre_Comun ,weight = 3, radius=4, 
              color=cf2(wildcats$Nombre_Comun), stroke = TRUE, fillOpacity = 0.5)  %>%
  addLegend("bottomright", colors= c("#ffa500", "#13ED3F"), labels=c("Puma", "Jaguar"), title="Especies") %>% 
  addLayersControl(
    baseGroups = c("OSM (default)", "Toner", "Toner Lite"),
    overlayGroups = wildcats$Nombre_Comun,
    options = layersControlOptions(collapsed = TRUE)
  )
m2

6.3.2 Selección final de la base de datos

6.3.2.1 Inclusion de nuevas variables

wildcats$Hora_Registro <- strftime(wildcats$Fecha_Registro, format="%H") #col7
wildcats$Fecha_Registro <- as.POSIXct(wildcats$Fecha_Registro, tz="", format="%Y-%m-%d") #col5
wildcats$Y.M.R <- strftime(wildcats$Fecha_Registro, format="%Y-%m") #col8
wildcats$M.R <- format(wildcats$Fecha_Registro, "%m") #col9
wildcats$Y.R <- format(wildcats$Fecha_Registro, "%Y") #col10

wildcats <- wildcats[, colnames(wildcats)[c(1:4,6,5,8,10,9,7)]] #Reorganizacion de las columnas

6.3.2.2 Visualizacion de la base de datos

datatable(wildcats, extensions = c('Buttons','FixedColumns', 'RowReorder'),
  filter = list(position = 'top', clear = FALSE),
  options = list(dom = 'Bfrtip',
    buttons = c('copy', 'csv', 'excel', 'pdf', 'print'), dom = 't',
    scrollX = TRUE,
    fixedColumns = TRUE,
    rowReorder = TRUE, order = list(c(0 , 'asc'))))

7 Resultados

7.1 Estadísticas espacio-temporales sobre los pumas y jaguares de un sector del continente americano

7.1.1 Especie por pais y mes

plot_ly(x= wildcats$Pais,split=wildcats$Nombre_Comun, color= I(cf2(wildcats$Nombre_Comun)), 
        frame=wildcats$M.R, type = "histogram", alpha = 0.6)%>%
  layout(barmode = "overlay")

7.1.2 Cantidad de Jaguares por pais

c.jaguars <- wildcats%>%
              dplyr::filter(Nombre_Comun == "Jaguar")%>%
              dplyr::group_by(Pais)%>%
              dplyr::summarise(Cant=n())
c.pumas <- wildcats%>%
              dplyr::filter(Nombre_Comun == "Puma")%>%
              dplyr::group_by(Pais)%>%
              dplyr::summarise(Cant=n())
data(worldgeojson, package = "highcharter")

highchart() %>% 
  hc_title(text = "Cantidad de Jaguares") %>% 
  hc_subtitle(text = "por pais") %>% 
  hc_add_series_map(map = worldgeojson, df = c.jaguars, name = "cantidad de individuos",
                    value = "Cant", joinBy = c("name", "Pais"), allAreas=FALSE, animation=TRUE)%>% 
  hc_colorAxis(min= 1,
            type='logarithmic',
            stops = color_stops())%>% 
  hc_tooltip(useHTML = TRUE, headerFormat = "",
             pointFormat = "Pais: {point.name}/ Cantidad de especies: {point.Cant}")%>%
  hc_legend(layout= 'horizontal',
            verticalAlign= 'top')%>%
  hc_mapNavigation(enabled = TRUE)

7.1.3 Cantidad de Pumas por pais

data(worldgeojson, package = "highcharter")

highchart() %>% 
  hc_title(text = "Cantidad de Pumas") %>% 
  hc_subtitle(text = "por pais") %>% 
  hc_add_series_map(map = worldgeojson, df = c.pumas, name = "cantidad de individuos",
                    value = "Cant", joinBy = c("name", "Pais"), allAreas= FALSE, animation=TRUE)%>% 
  hc_colorAxis(min= 1,
            type='logarithmic',
            stops = color_stops())%>% 
  hc_tooltip(useHTML = TRUE, headerFormat = "",
             pointFormat = "Pais: {point.name}/ Cantidad de especies: {point.Cant}")%>%
  hc_legend(layout= 'horizontal',
            verticalAlign= 'top')%>%
  hc_mapNavigation(enabled = TRUE)

7.2 Distribucion de las especies y lineas de tiempo

7.2.1 Especies por mes y año

ct.wildcats <- wildcats%>%
              dplyr::group_by(Fecha_Registro, Nombre_Comun)%>%
              dplyr::summarise(Cant=n())

ct.wildcats <-reshape(ct.wildcats, idvar = "Fecha_Registro", timevar="Nombre_Comun", direction = "wide")
highchart(type = "stock") %>% 
  hc_title(text = "Cantidad de especies") %>% 
  hc_subtitle(text = "en la linea de tiempo de observacion") %>% 
  hc_add_series(wildcats$Y.M.R, type = "spline", pointInterval= 30*24*3600*1000, name="Cant de especies en la linea de tiempo")%>%
  hc_xAxis(title= list(text='Fechas'), type='datetime')%>%
  hc_yAxis(title= list(text='Cantidad de especies observadas'), opposite = FALSE)

7.2.2 Distribucion de las especies por mes

e.m <- ggplot(wildcats, aes(x=Longitud, y=Latitud, frame=M.R, text = paste("Pais:", Pais))) +
  borders("world", xlim = c(-120, -30), ylim = c(-50, 25))+
  geom_point(data=filter(wildcats), shape = 21, colour = cf2(wildcats$Nombre_Comun), size = 1)+
  facet_grid(cols = vars(Nombre_Comun)) +
  labs(title = "Distribucion de las especies observadas segun mes")

ggplotly(e.m)%>% 
  animation_opts(
    2000, easing = "elastic", redraw = FALSE
  )%>%
  animation_slider(
    currentvalue = list(prefix = "Mes ", font = list(color="red"))
  )

7.2.3 Patron geometrico de distribucion

pd <- qplot(data = wildcats, x=Longitud, y=Latitud) +
  stat_ellipse(geom = "polygon", alpha = 1/2, aes(fill = Nombre_Comun, frame=M.R))+
  borders("world", xlim = c(-120, -30), ylim = c(-50, 25))+
  geom_point(data=filter(wildcats), shape = 21, colour = cf2(wildcats$Nombre_Comun), size = 0.5, alpha=0.3)+
  facet_grid(cols = vars(Nombre_Comun)) +
  labs(title = "Patron de Distribucion de las Especies")

ggplotly(pd)

7.3 Analisis de densidad espacial de las especies

7.3.1 Densidad de cada una de las especies

dn <- ggplot(wildcats, aes(x=Longitud, y=Latitud)) +
  stat_density2d(aes(fill = stat(level)), geom="polygon") +
  borders("world", regions=".", xlim = c(-120, -30), ylim = c(-50, 25))+
  geom_point(data=filter(wildcats), shape = 21, colour = cf2(wildcats$Nombre_Comun), size = 0.5, alpha=0.3)+
  facet_grid(cols = vars(Nombre_Comun)) +
  scale_fill_viridis_c(option = "viridis") +
  theme(legend.position = "magma") +
  labs(title = "Densidad segun distribucion por especie")

ggplotly(dn)

7.3.2 Interpolacion